NASA/TM— 2009-215301 


AIAA-2008^1743 



Fan Flutter Computations Using 
the Harmonic Balance Method 


MilindA. Bakhle 

Glenn Research Center, Cleveland, Ohio 
Jeffrey P. Thomas 

Duke University, Durham, North Carolina 
T.S.R. Reddy 

University of Toledo, Toledo, Ohio 


February 2009 


NASA STI Program ... in Profile 


Since its founding, NASA has been dedicated to the 
advancement of aeronautics and space science. The 
NASA Scientific and Technical Information (STI) 
program plays a key part in helping NASA maintain 
this important role. 

The NASA STI Program operates under the auspices 
of the Agency Chief Information Officer. It collects, 
organizes, provides for archiving, and disseminates 
NASA’s STI. The NASA STI program provides access 
to the NASA Aeronautics and Space Database and 
its public interface, the NASA Technical Reports 
Server, thus providing one of the largest collections 
of aeronautical and space science STI in the world. 
Results are published in both non-NASA channels 
and by NASA in the NASA STI Report Series, which 
includes the following report types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant phase 
of research that present the results of NASA 
programs and include extensive data or theoretical 
analysis. Includes compilations of significant 
scientific and technical data and information 
deemed to be of continuing reference value. 
NASA counterpart of peer-reviewed formal 
professional papers but has less stringent 
limitations on manuscript length and extent of 
graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 
and technical findings that are preliminary or 
of specialized interest, e.g., quick release 
reports, working papers, and bibliographies that 
contain minimal annotation. Does not contain 
extensive analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


papers from scientific and technical 
conferences, symposia, seminars, or other 
meetings sponsored or cosponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from 
NASA programs, projects, and missions, often 
concerned with subjects having substantial 
public interest. 

• TECHNICAL TRANSLATION. English- 
language translations of foreign scientific and 
technical material pertinent to NASA’s mission. 

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results. 

For more information about the NASA STI 

program, see the following: 

• Access the NASA STI program home page at 
http://www.sti. nasa.gov 

• E-mail your question via the Internet to help@ 
sti.nasa.gov 

• Fax your question to the NASA STI Help Desk 
at 301-621-0134 

• Telephone the NASA STI Help Desk at 
301-621-0390 

• Write to: 

NASA Center for AeroSpace Information (CASI) 
7115 Standard Drive 
Hanover, MD 21076-1320 


CONFERENCE PUBLICATION. Collected 


NASA/TM— 2009-215301 


AIAA-2008^1743 



Fan Flutter Computations Using 
the Harmonic Balance Method 


MilindA. Bakhle 

Glenn Research Center, Cleveland, Ohio 
Jeffrey P. Thomas 

Duke University, Durham, North Carolina 
T.S.R. Reddy 

University of Toledo, Toledo, Ohio 


Prepared for the 

44th Joint Propulsion Conference and Exhibit 
cosponsored by AIAA, ASME, SAE, and ASEE 
Hartford, Connecticut, July 21-23, 2008 


National Aeronautics and 
Space Administration 


Glenn Research Center 
Cleveland, Ohio 44135 


February 2009 


Acknowledgments 


Support for this work, from the NASA Fundamental Aeronautics Program, Subsonic Fixed Wing Project 
and Supersonics Project, is greatfully acknowledged. 


This work was sponsored by the Fundamental Aeronautics Program 
at the NASA Glenn Research Center. 


Level of Review : This material has been technically reviewed by technical management. 


Available from 


NASA Center for Aerospace Information 
7115 Standard Drive 
Hanover, MD 21076-1320 


National Technical Information Service 
5285 Port Royal Road 
Springfield, VA 22161 


Available electronically at http://gltrs.grc.nasa.gov 


Fan Flutter Computations Using the Harmonic Balance Method 


Milind A. Bakhle 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 

Jeffrey P. Thomas 
Duke University 
Durham, North Carolina 27708 

T. S. R. Reddy 
University of Toledo 
Toledo, Ohio 43606 


Abstract 

An experimental forward-swept fan encountered flutter at part-speed conditions during wind tunnel testing. A 
new propulsion aeroelasticity code, based on a computational fluid dynamics (CFD) approach, was used to model 
the aero elastic behavior of this fan. This three-dimensional code models the unsteady flowfield due to blade 
vibrations using a harmonic balance method to solve the Navier- Stokes equations. This paper describes the flutter 
calculations and compares the results to experimental measurements and previous results from a time- accurate 
propulsion aeroelasticity code. 


I. Introduction 

Research has been on-going in the development, validation and application of high-fidelity models for 
aeroelastic vibrations in aircraft engine fan, compressor, and turbine blades. Recent work has included time-domain 
solution of the Reynolds- averaged Navier- Stokes (RANS) equations to provide the unsteady flowfield and unsteady 
aerodynamic forces on the blades. An example of such work is the TURBO aeroelastic analysis code (refs. 1 and 2). 
Such high-fidelity time-domain models require large numbers of computations and take a long time for startup 
transients to decay before the final periodic solution is obtained. An alternate approach (ref. 3) is to use the 
periodicity in time of typical turbomachinery flows to represent each flow variable by a Fourier series in time, 
leading to a harmonic balance form of the Navier-Stokes equations. Solutions to these equations can be obtained 
using methods that are typically used for steady flow problems such as pseudo-time marching and local time 
stepping. Thus, the harmonic balance approach leads to a method that is significantly faster than the typical time- 
domain solution method (ref. 3). 

In the present study, the harmonic balance aeroelastic analysis method and computer code developed by Hall et 
al. (refs. 3 and 4) is applied to the problem of part-speed fan flutter. The configuration selected is an experimental 
fan for which wind-tunnel data is available, and for which aeroelastic computations have been previously performed 
using a time-domain RANS aeroelastic code (ref. 1). The design and testing of the advanced single-stage transonic 
fan has been described in detail in a recent report (ref. 5). This fan was designed with aggressive goals for 
performance and noise reduction. During rig testing in a high-speed wind tunnel, the fan performed well at design 
speed, and was successfully throttled to the stall line. However, flutter was encountered just above the operating line 
at part-speed conditions. The flutter mode was identified as the first bending mode of the airfoil, in a two nodal 
diameter forward-traveling wave pattern. 

In the present work, a computational study was performed to examine the flutter characteristics using a new 
state-of-the-art propulsion aeroelasticity code based on the harmonic balance method. This code (ref. 4) solves the 
three-dimensional unsteady, Reynolds- Averaged Navier-Stokes equations with the ability to model a rotating blade 
row with harmonic blade vibrations or incoming periodic distortions. For flutter calculations, the blade vibration is 
prescribed to be the modal deflection and a frequency, both of which are calculated from a separate structural 
dynamics analysis. Computations are performed in a single blade passage for both steady and unsteady analysis. The 
resulting unsteady pressures are used to calculate a generalized aerodynamic force and then an eigenvalue problem 
is solved to calculate the aerodynamic damping, which is used to determine stability of the blade in a specific 
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vibration mode. Note that in the present study, the harmonic balance aeroelastic code was run with very small 
amplitude vibrations of the fan blade to calculate linearized unsteady aerodynamics for a conventional linear flutter 
analysis. Although the harmonic balance method provides a computationally efficient frequency-domain solution for 
modeling nonlinear unsteady aerodynamic effects due to finite amplitude blade vibrations motions, the present 
calculations have been restricted to small amplitude vibrations in order to examine only linearized aerodynamic 
characteristics. 


II. Analysis 

A. Aeroelastic Model 

The equations of motion for a fan blade (with all blades assumed to be identical) can be written as 

(i) 

where [M] and [X] are generalized mass and stiffness matrices, {q} is the generalized displacement vector, and [A] 
is the blade vibration-dependent generalized aerodynamic force matrix. The matrices [ M ], [X] and [A] are of size 
NMxNM; {g}is of size NMx 1; NM is the number of modes. 

The elements of [M] and [K\ are obtained from a free-vibration analysis using a commercial structural dynamics 
analysis software. The matrices [M] and [X] are diagonal and their non-zero elements are related as 

K , =M l m i 2 (] + 2iQ (2) 

where (0/ is the natural frequency of the i th mode, and C, is the structural damping ratio; usually the mode shapes are 
mass-normalized and therefore M t = 1 . 

Since all the blades are identical (that is, a tuned rotor), the aeroelastic modes consist of individual blades 
vibrating with equal amplitudes at a fixed interblade phase angle between adjacent blades. Hence, the motion of the 
s th blade in r th interblade phase angle mode can be written as 


{q s } = {q r }e i ™ e*'* (3) 

where co is the vibration frequency, o r is the interblade phase angle related to nodal diameter (ND) pattern of the 
traveling wave and number of blades Abides as 


O r 2lZ AD/Xbiades 


( 4 ) 


Thus, the equations of motion for a blade become 

-CO 2 [M]{q r }+[K] {q r }=\A r ]{q r } (5) 

B. Calculation of Aerodynamic Force Matrix [A r \ 

Flutter analysis requires the calculation of all elements of the generalized aerodynamic force matrix [ A r ]. 
Unsteady flowfield computations are carried out for each vibration mode and assumed frequency. For a selected 
value of interblade phase angle, the harmonic balance code is used to calculate the unsteady pressure distribution on 
the blade surface, which is further used to calculate the (complex-valued) elements of the generalized aerodynamic 
force matrix [A r \. This calculation is repeated for AWdes interblade phase angles given by eq. (4). 

C. Flutter Analysis 

An eigenvalue approach is used to determine flutter stability. Equation (5) is written in a standard eigenvalue 
form as: 
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[M-y[e]]K}={o} 


( 6 ) 


where [P] =([A^]-[A r ])/cOp; [2] = [M]; y = (co/co 0 ) 2 ; and co 0 is the assumed frequency used in the calculation of 
the elements of the aerodynamic force matrix [A r ]. 

The solution of the above eigenvalue problem results in NM complex eigenvalues of the form 

i(m/m 0 )=i-y[y = ^±iv (7) 

The real part of the eigenvalue ( ft ) represents the damping ratio, and the imaginary part ( v ) represents the 
damped frequency; flutter occurs if jl > 0 for any eigenvalue. 

For turbomachinery blades, flutter analysis is usually carried out with a single vibration mode. In this case, 
eq. (6) reduces to a scalar equation and the eigenvalues can be written as 

Y = P/Q (8) 

In the present work, the structural damping is set to zero and the vibration mode is mass-normalized, leading to 
the following simplified solution for the eigenvalue problem 

Y = (® natural ~ ^)/®o ( 9 ) 

Note that the frequency used in the unsteady aerodynamic calculations (co 0 ) has the same value as the blade 
natural frequency (co na turai)- Also note that since the structural damping is set to zero, the damping is referred to as 
aerodynamic damping in the following section, but it is opposite in sign to (I . 


III. Results 

In this section, the results of the steady and unsteady computations are presented. Steady and unsteady 
computations were carried out using the TURBO code at rotational speeds of 85 and 75 percent. These two speeds 
were selected since during testing, large vibratory responses due to flutter were encountered on the stall side of the 
operating line at these speeds. Note that during testing at 100 percent (design) speed, no flutter was encountered. 

Note that in the present work, the harmonic balance aeroelastic code was run with very small amplitude 
vibrations of the fan blade to calculate linearized unsteady aerodynamics for a conventional linear flutter analysis. 

A. Steady Computational Results 

In this study, the airfoil geometry used was calculated based on static deflections obtained from structural 
analysis for 85 percent speed. The static deflections included the effects of rotational speed, applied pressures, and 
blade temperatures; nonlinear geometry effects were also included in the analyses. The blade geometry at 85 percent 
speed and nominal operating conditions was used for all computations. Previous calculations (ref. 1) for this case 
have shown that changes from the nominal blade geometry due to changes in rotational speed were not significant 
and therefore a single geometry was used for computations at all speeds and operating conditions in the present 
work. The computational grid used was generated using commercial software. The grid is shown in figure 1 ; the grid 
size is 193x33x49 for the block that wraps around the blade airfoil with 193 grid points around the airfoil, 33 grid 
points in the circumferential direction, and 49 grid points in the spanwise direction. The grid blocks in the inlet and 
exit sections are each 17x33x49 with 17 points in the streamwise direction, 33 grid points in the circumferential 
direction, and 49 grid points in the spanwise direction. The tip clearance was based on rig test measurements and is 
modeled using 9 points between the blade tip and casing. 
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Figure 1. — Computational mesh. 




Figure 2. — Steady solution convergence. 


Figure 3. — Fan map showing current 
performance predictions. 


The inlet flow conditions used in the computations consisted of circumferentially-averaged radial profiles of 
total pressure, total temperature, and flow angles. These prescribed profiles were based on rig measurements, 
supplemented by previous steady computations. The exit flow conditions consisted of a circumferentially-averaged 
radial profile of static pressure. This profile, which was based on pressure measurements at design speed, was used 
with uniform scaling for the computations at all speeds and at all conditions. 

Figure 2 shows the convergence of the steady computations as a plot of non-dimensional torque with iteration 
counter (100 iterations per counter). The computations were carried out at 85 percent speed for various values of 
imposed back-pressure at the exit boundary. As can be noted from figure 2, excellent convergence was obtained for 
all operating conditions. Similar convergence was obtained for the computations at 75 percent speed. All steady 
results presented here are from well-converged solutions. 

The experimental fan map is shown in figure 3 along with the results of current steady computations denoted as 
HB. Several computations were done for 85 and 75 percent speeds and the results correlate well with measurements. 
Of particular interest is the last condition computed on the stall side, beyond which no converged solutions were 
obtained. The mass flow rate at the near-stall point is in excellent agreement with measurements; the predicted 
pressure ratio is approximately 3 percent lower at both speeds. No converged solutions were obtained at 85 percent 
speed for lower mass flow rates beyond those shown in figure 3. For 75 percent speed, additional calculations are 
required with small increases in back-pressure to determine if additional stable solutions can be obtained for lower 
mass flow rates. Note that the stall line shown in figure 3 was estimated based on prior experience and was not 
confirmed during rig testing because of the occurrence of flutter. Also, since the stall side was of primary interest, 
no attempt was made to determine the upper bound on choke flow. 
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Inlet Corrected Flow 


Figure 4. — Fan map showing current performance 
predictions compared with TURBO results (ref. 1). 


Figure 4 shows a comparison of the present results for performance compared to the previous computational 
results from the TURBO code (ref. 1). The present results show good correlation with the TURBO results, except 
that both 75 and 85 percent speed lines are shifted towards lower mass flow rate and lower pressure ratio. The 
present calculations provide converged results for lower mass flow rates than were obtained with TURBO. Also, the 
present results consistently show a lower pressure ratio as compared to TURBO and the test data. These differences 
need to be investigated further. It should be noted that significant differences exist between the numerical method 
used in the previous TURBO computations and the current computations with the harmonic balance code. These 
differences include time-domain versus frequency-domain, numerical discretization, grid, and algorithm used to 
solve the RANS equations, turbulence modeling, and others. Therefore, the differences in results noted in figure 3 
are not entirely surprising. Also, additional sensitivity studies are required to determine the variability in the newly 
computed results. 


B. Unsteady Computational Results 

Unsteady computations were performed using the harmonic balance code at 85 and 75 percent speeds to predict 
the aeroelastic stability of the blade. Several nodal diameters covering the entire possible range were considered. 
The steady condition near peak-efficiency was used along with a condition near the stall line. This allowed the trend 
along the speed line to be assessed. Note that the results presented in this paper were all computed with only the 0 th 
and 1 st harmonics included in the aerodynamic analysis using the harmonic balance method. 

Figure 5 shows the convergence of the unsteady computations as a plot of non-dimensional generalized force 
(real and imaginary parts) with iteration counter (100 iterations per counter) for various nodal diameters. These 
computations were for 85 percent speed and for an operating point near the peak-efficiency condition with a non- 
dimensional back-pressure of 1.0. As can be noted, although the rates of convergence vary with nodal diameter, all 
the results are well-converged. Similar convergence characteristics were observed at the near-stall condition with a 
non-dimensional back-pressure of 1 .05. 

A sensitivity study was performed to investigate the effect of the unsteady grid scaling parameter. The nominal 
value of this numerical input parameter was increased by factors of 2 and 5, and then decreased by the same factors. 
There was no change in the generalized force results within plotting accuracy — demonstrating that for the present 
calculations, the prescribed unsteady grid scaling has no effect on the results. Recall that in the present study, the 
harmonic balance aeroelastic code was run with very small amplitude vibrations of the fan blade to calculate 
linearized unsteady aerodynamics for a conventional linear flutter analysis. 

Figure 6 shows the variation of the converged generalized force with nodal diameter for the two different values 
of back-pressure at 85 percent speed. Note the significant variation in the imaginary part of the generalized force. 
For a single vibration mode, the stability is determined by the imaginary part of the generalized force and figure 6 
shows that the stability varies significantly with nodal diameter. Further, with the change in back-pressure from 
peak-efficiency towards stall, it can be seen that the imaginary part of the generalized force drops closer to zero at 


NASA/TM — 2009-2 15301 


5 


ND= 2. This trend is seen more clearly in figure 7. The imaginary part of the generalized force is seen to clearly 
move towards zero as the back-pressure increases — going from peak-efficiency towards stall. Note that figure 7 
includes results at non-dimensional back-pressure values of 1.052 and 1.0522, with a change in the sign of the 
imaginary part occurring between these conditions. 




Figure 5. — Unsteady solution convergence. Figure 6. — Generalized force variation 

with nodal diameter. 



Non-Dimensional Back-Pressure 

Figure 7. — Variation of generalized force (imaginary part) 
with back-pressure for ND=2. 
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Figure 8. — Variation of aerodynamic damping with 
nodal diameter at 85 percent speed. 




Non-Dimensional Mass Flow Rate Non-Dimensional Mass Flow Rate 


Figure 9. — Variation of aerodynamic damping with 
non-dimensional mass flow rate for ND=2 at 
85 percent speed. 


Figure 10. — Variation of aerodynamic damping with 
non-dimensional mass flow rate for ND=2 
at 75 percent speed. 


The generalized force was used to calculate the aerodynamic damping, which is plotted in figure 8. The results 
are presented for back-pressure values near peak-efficiency and near-stall conditions. The aerodynamic damping is 
seen to vary by an order of magnitude with nodal diameter of the traveling wave, with the minimum occurring at 
ND= 2 (forward traveling wave). Note that although the damping value is nearly zero at the NS condition, it is still 
positive even at its minimum value, indicating stability at a non-dimensional back-pressure of 1 .05. 

Additional calculations were carried out for a small increase in back-pressure to identify the condition of zero 
aerodynamic damping (flutter). In the absence of these additional results at back-pressure of 1.052 and 1.0522, it 
would be necessary to extrapolate the damping as a function of mass flow to identify the flutter point. Figure 9 
shows the variation of aerodynamic damping with mass flow rate with a clear monotonic drop in stability along the 
speed line towards stall. The results of unsteady calculations at 75 percent speed show trends that are nearly the 
same as those in figures 5 to 9. Figure 10 shows the corresponding variation of aerodynamic damping along the 
75 percent speed line. Note that the aerodynamic damping is clearly negative at the lowest mass flow rate condition 
and the flutter point is obtained by interpolation. 
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Inlet Corrected Flow 

Figure 1 1 . — Fan map showing calculated flutter point 
along with the measured flutter boundary. 


The calculated flutter points at 85 and 75 percent speeds are plotted on the fan map along with the measured 
flutter boundary in figure 1 1 . It can be seen that the calculated flutter point is very close to the measured value at 
85 percent sped with a difference of less than 0.5 percent in mass flow rate. As noted previously, the pressure ratio is 
under-predicted and the difference from the measured data is 3 percent at the flutter point. For comparison, the 
flutter point calculated using TURBO is also plotted. Both calculations are very close to the test data; the present 
calculated result is slightly closer to the measurement in mass flow rate, but farther in pressure ratio as compared to 
the TURBO result. 

At 75 percent speed, the difference between the present results and measurements is higher (about 4 percent) 
than at 85 percent speed. The TURBO result also shows a larger difference from measurement at 75 percent speed 
than at 85 percent speed. 


IV. Conclusion 

Steady and unsteady computations have been carried out using a Propulsion Aeroelasticity code based on the 
harmonic balance method of solving the Reynolds- Averaged Navier-Stokes equations. The configuration selected 
was an experimental fan that encountered flutter during wind tunnel testing. The computational results were 
summarized and compared with experimental data and results from a previous computational study with a time- 
domain propulsion aeroelasticity code (TURBO). Overall, the results are in good agreement with experimental data 
and the previous TURBO results. The steady results were compared on the performance map and showed good 
correlation with data and TURBO results. However, the present results under-predict the pressure ratio. Also, the 
present analysis provides solutions at slightly lower mass flow rates than were obtained with TURBO. The 
systematic under-prediction of the pressure ratio by about 3 percent needs to be investigated further. A detailed look 
at the flowfield will provide additional information regarding the source of the lower pressure ratio. Also, numerical 
studies are required to investigate the sensitivity of the present results to numerical input parameters. 

The flutter results correlated very well with the experimental data and the previous results from TURBO. The 
least stable nodal diameter pattern (interblade phase angle) matched correctly with the experimental observations 
and TURBO results. The calculated flutter condition agreed closely with the experimental data and TURBO results. 
A sensitivity study on the grid motion scaling parameter showed that the present results were insensitive to 
variations in the selected value of this parameter, thus demonstrating the linearity of the results for the selected small 
amplitude of blade vibration. 

In the present work, computations were restricted to the first bending mode of vibration. Computations with 
other modes are required to show that the first mode is the least stable mode. Additional steady and unsteady 
computations are required to determine at 100 percent speed to demonstrate that no flutter is predicted for 
100 percent speed, consistent with experimental observations and previous calculations with TURBO. 

The real power of the harmonic balance method is that it can be a computationally efficient frequency- domain 
alternative to more computationally expensive time-domain methods for modeling nonlinear unsteady aerodynamic 
effects brought about by finite- amplitude blade vibrations. However, the present nonlinear harmonic balance 
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aeroelastic code also provides capabilities for a linearized analysis as a special case. In the present study, very small 
amplitude vibrations were prescribed to calculate linearized unsteady aerodynamics for a conventional linear flutter 
analysis as a first step to understand the characteristics of the harmonic balance code and of the flutter of an 
experimental fan. Future work will investigate the non-linear amplitude- dependent effects using the same harmonic 
balance aeroelastic code. 
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